#include "yocto/swamp/common.hpp"
#include "yocto/string/vfs-utils.hpp"
#include "yocto/lua/lua-state.hpp"
#include "yocto/lua/lua-config.hpp"
#include "yocto/ios/ocstream.hpp"
#include "yocto/code/utils.hpp"
#include "yocto/math/types.hpp"

#include "yocto/math/ode/stiff-drvkr.hpp"

#include "yocto/wtime.hpp"
#include "yocto/string/conv.hpp"
#include "yocto/auto-ptr.hpp"

#include "yocto/z++/gzstream.hpp"
#include "yocto/eta.hpp"
#include "yocto/duration.hpp"
#include "yocto/string/env.hpp"
#include "yocto/threading/thread.hpp"
#include "yocto/hw.hpp"

using namespace yocto;
using namespace swamp;
using namespace math;
using namespace threading;

static lua_State *L = NULL;

typedef double                         Real;
typedef layout1D                       Layout;
typedef coord1D                        Coord;
typedef fields_setup<Layout>           FieldsSetup;
typedef array1D<Real>                  Array1D;
typedef ghosts_setup<Coord>            GhostsSetup;
typedef workspace<Layout, Real, rmesh> Workspace;
typedef region1D<Real>::type           RegionType;
typedef vertex1D<Real>::type           Vertex;
typedef ode::stiff_drvkr<Real>::type   ImplicitIntegrator;

typedef ode::field<Real>::type         DiffEq;
typedef ode::field<Real>::diff         Jacobn;

static const Real ALPHA_SAFE = 0.1;
static const Real dt0        = 1e-6;


static bool GFX = true;
static void gfxDraw();


////////////////////////////////////////////////////////////////////////////////
//
// Fields and Ghosts Setup for the simulation
//
////////////////////////////////////////////////////////////////////////////////
#define SIM_CFG(NAME) NAME( Lua::Config::Get<lua_Number>(L,#NAME) )

class SimParam : public FieldsSetup
{
public:
    const GhostsSetup sim_ghosts;
    const size_t      volumes;
    const Layout      sim_layout;
    const Real        length;
    const Real        Dh;
    const Real        Dw;
    const Real        Kw;
    Real              tol;
    const Real        pH_left;
    const Real        pH_right;
    const Real        kd;
    const Real        k2;
    
    explicit SimParam( lua_State *L) : 
    sim_ghosts(),
    SIM_CFG(volumes),
    sim_layout(0,volumes+1),
    SIM_CFG(length),
    SIM_CFG(Dh),
    SIM_CFG(Dw),
    SIM_CFG(Kw),
    tol( 1e-3 ),
    SIM_CFG(pH_left),
    SIM_CFG(pH_right),
    SIM_CFG(kd),
    SIM_CFG(k2)
    {
        Y_SWAMP_DECL_VAR( *this, "h",    Array1D);
        Y_SWAMP_DECL_VAR( *this, "w",    Array1D);
        Y_SWAMP_DECL_VAR( *this, "Fh",   Array1D);
        Y_SWAMP_DECL_VAR( *this, "Fw",   Array1D);
        Y_SWAMP_DECL_VAR( *this, "ih",   Array1D);
        Y_SWAMP_DECL_VAR( *this, "iw",   Array1D);
        Y_SWAMP_DECL_VAR( *this, "step", Array1D);
    }
    
    
    virtual ~SimParam() throw()
    {
    }
    
    
private:
    YOCTO_DISABLE_COPY_AND_ASSIGN(SimParam);
};


////////////////////////////////////////////////////////////////////////////////
//
// Simulation with its data
//
////////////////////////////////////////////////////////////////////////////////
class Simulation : public SimParam, public Workspace
{
public:
    const size_t       N;
    const Real         dx;
    const Real         inv_dx;
    const Array1D     &X;
    Array1D           &h;
    Array1D           &w;
    Array1D           &Fh;
    Array1D           &Fw;
    Array1D           &ih;
    Array1D           &iw;
    Array1D           &step;
    ImplicitIntegrator imODE;
    DiffEq             drvs;
    Jacobn             djac;
    vector<Real>       Y;
    size_t             num_newton;
    Real               ave_newton;
    wtime              chrono;
    double             t_diff;
    double             t_proj;
    
#define SIM_LINK(NAME) NAME( (*this)[ #NAME ].as<Array1D>() )
    
    explicit Simulation( lua_State *L ) : 
    SimParam(L),
    Workspace( sim_layout, sim_ghosts, *this ),
    N( upper ),
    dx( length/N ),
    inv_dx( N/length ),
    X( mesh.X() ),
    SIM_LINK(h),
    SIM_LINK(w),
    SIM_LINK(Fh),
    SIM_LINK(Fw),
    SIM_LINK(ih),
    SIM_LINK(iw),
    SIM_LINK(step),
    imODE(),
    drvs( this, & Simulation:: diffeq ),
    djac( this, & Simulation:: jacobn ),
    Y(2,0.0),
    num_newton(0),
    ave_newton(0),
    t_diff(0),
    t_proj(0)
    {
        //----------------------------------------------------------------------
        // create mesh
        //----------------------------------------------------------------------
        const RegionType r(0,length);
        mesh.regular_map_to(r,*this);
        
        //----------------------------------------------------------------------
        // prepare solver
        //----------------------------------------------------------------------
        imODE.start(2);
        
        chrono.start();
    }
    
    virtual ~Simulation() throw()
    {
    }
    
    ////////////////////////////////////////////////////////////////////////////
    //
    //! initialize with boundary pH
    //
    ////////////////////////////////////////////////////////////////////////////
    void initialize() throw()
    {
        h[0]    = pow(10,-pH_left);
        w[0]    = Kw/h[0];
        step[0] = dt0;
        const Real h_right = pow(10,-pH_right);
        const Real w_right = Kw/h_right;
        for( size_t i=1; i <=N; ++i )
        {
            h[i]    = h_right;
            w[i]    = w_right;
            step[i] = dt0;
        }
        //h[N] = h[0]/2; w[N] = Kw/h[N];
        t_diff     = 0;
        t_proj     = 0;
        num_newton = 0;
        ave_newton = 0;
        imODE.eps  = tol;
    }
    
    ////////////////////////////////////////////////////////////////////////////
    //
    //! compute flux for one specie
    //
    ////////////////////////////////////////////////////////////////////////////
    void compute_flux( Array1D &F, const Array1D &C, const Real D ) const
    {
        const Real fac  = D * inv_dx;
        for( size_t i=0; i < N; ++i )
        {
            F[i] = - fac * (C[i+1] - C[i]);
        }        
    }
    
    ////////////////////////////////////////////////////////////////////////////
    //
    //! compute all fluxes
    //
    ////////////////////////////////////////////////////////////////////////////
    void compute_fluxes( )
    {
        compute_flux( Fh, h, Dh );
        compute_flux( Fw, w, Dw );
    }
    
    ////////////////////////////////////////////////////////////////////////////
    //
    //! compute increase for one species
    //
    ////////////////////////////////////////////////////////////////////////////
    void compute_increase( const Real dt, Array1D &inc, const Array1D &F ) const
    {
        const Real fac = dt * inv_dx;
        for(size_t i=1; i < N; ++i )
        {
            inc[i] = -fac*( F[i]-F[i-1] );
        }
    }
    
    ////////////////////////////////////////////////////////////////////////////
    //
    //! compute all increases
    //
    ////////////////////////////////////////////////////////////////////////////
    void compute_increases( const Real dt )
    {
        compute_increase( dt, ih, Fh );
        compute_increase( dt, iw, Fw );
    }
    
    ////////////////////////////////////////////////////////////////////////////
    //
    //! compute full diffusive increase
    //
    ////////////////////////////////////////////////////////////////////////////
    void plain_diffusion( const Real dt )
    {
        compute_fluxes();
        compute_increases(dt);
    }
    
    
    
    ////////////////////////////////////////////////////////////////////////////
    //
    //! plain update with implicit reaction
    //
    ////////////////////////////////////////////////////////////////////////////
    void implicit_plain_update( const Real dt)
    {
        for( size_t i=1; i < N; ++i )
        {
            Y[1] = ( h[i] += ih[i]);
            Y[2] = ( w[i] += iw[i]);
            imODE( drvs, djac, Y, 0, dt, step[i] );
            h[i] = Y[1];
            w[i] = Y[2];
        }
    }
    
    ////////////////////////////////////////////////////////////////////////////
    //
    //! compute relaxed diffusion
    //
    ////////////////////////////////////////////////////////////////////////////
    void relaxed_diffusion( const Real dt )
    {
        plain_diffusion(dt);
        for( size_t i=1; i < N; ++i )
        {
            const Real H  = h[i];
            const Real W  = w[i];
            const Real dH = ih[i];
            const Real dW = iw[i];
            const Real fac = 1.0/(H+W);
            ih[i] = H * ( dH - dW ) * fac;
            iw[i] = W * ( dW - dH ) * fac;
        }
    }
    
    ////////////////////////////////////////////////////////////////////////////
    //
    //! apply increases and perform newton step
    //
    ////////////////////////////////////////////////////////////////////////////
    void relaxed_update(const Real)
    {
        num_newton = 0;
        for( size_t i=1; i < N; ++i )
        {
            h[i] += ih[i];
            w[i] += iw[i];
            newton( h[i], w[i] );
        }
        ave_newton = num_newton / Real(N-1);
    }
    
    //! local newton step to the provided tolerance
    void newton( Real & H, Real &W )
    {
        assert(H>0);
        assert(W>0);
        while( true )
        {
            ++num_newton;
            const Real Gam  = Kw - (H*W);
            const Real dxi  = Gam/(H+W);
            const Real dh   = dxi;
            const Real dw   = dxi;
            H += dh;
            W += dw;
            if( Fabs(dh) <= Fabs( tol * H ) && Fabs(dw) <= Fabs( tol * W ) )
                break;
        }
    }
    
    
    
    //! rescale increases
    inline void rescale_increases() throw()
    {
        for( size_t i=1; i <N; ++i )
        {
            ih[i] *= 0.5;
            iw[i] *= 0.5;
        }
    }
    
    
#define __CHRONO(tmx,CODE) do { const double t1 = chrono.query(); CODE; tmx += chrono.query() - t1; } while(false)
    //! adaptive diffusion 
    void adaptive_diffusion(Real dt, 
                            void (Simulation::*diffusion_step)(const Real),
                            void (Simulation::*update_step)(const Real)
                            )
    {
        do
        {
            //------------------------------------------------------------------
            //-- try full step 
            //------------------------------------------------------------------
            __CHRONO(t_diff,(this->*diffusion_step)(dt));
            
            //------------------------------------------------------------------
            //-- check what is possible
            //------------------------------------------------------------------
            Real dt_done = dt;
        CHECK_STEP:
            for( size_t i=1; i <= N; ++i )
            {
                assert( h[i]>0 );
                assert( w[i]>0 );
                if( ih[i] <= -h[i] || iw[i] <= -w[i] )
                {
                    dt_done *= 0.5;
                    rescale_increases();
                    goto CHECK_STEP;
                }
            }
            
            //------------------------------------------------------------------
            //-- update time left
            //------------------------------------------------------------------
            dt -= dt_done;
            
            //------------------------------------------------------------------
            //-- validate so far
            //------------------------------------------------------------------
            __CHRONO(t_proj,(this->*update_step)(dt_done));
        }
        while(dt>0);
        
    }
    
    
    
    void implicit_plain_scheme( Real dt )
    {
        adaptive_diffusion(dt, & Simulation::plain_diffusion, &Simulation::implicit_plain_update);
    }
    
    void relaxed_scheme( Real dt )
    {
        adaptive_diffusion(dt, & Simulation::relaxed_diffusion, &Simulation::relaxed_update );
    }
    
    
    void save_field( const char *filename, const Array1D &u, size_t nmin, size_t nmax )
    {
        assert(nmin<=nmax);
        assert(nmax<=N);
        ios::ocstream fp( filename, false );
        for( size_t i=nmin; i <= nmax; ++i )
        {
            fp("%g %g\n", X[i], u[i] );
        }
    }
    
    void save_perf( const string &filename, const Real t, bool append) const
    {
        ios::ocstream fp( filename, append );
        if(!append)
            fp("#t t_diff t_proj total newton\n");
        fp( "%g %g %g %g %g\n", t, t_diff, t_proj, t_diff+t_proj, ave_newton );
    }
    
    void save_db( const string &filename,const Real t,bool append) const
    {
        ios::ogzstream fp( filename, 9, append );
        
        //----------------------------------------------------------------------
        // save global info
        //----------------------------------------------------------------------
        fp.emit<double>(t);
        fp.emit<double>(t_diff);
        fp.emit<double>(t_proj);
        fp.emit<double>(ave_newton);
        
        //----------------------------------------------------------------------
        // save fields
        //----------------------------------------------------------------------
        fp.save( h.entry, h.bytes);
        fp.save( w.entry, w.bytes);
    }
    
    bool load_db( Real &t, ios::istream &fp )
    {
        char C = 0;
        if( ! fp.query(C) )
            return false;
        fp.store(C);
        //----------------------------------------------------------------------
        // load global info
        //----------------------------------------------------------------------
        t          = fp.read<double>();
        t_diff     = fp.read<double>();
        t_proj     = fp.read<double>();
        ave_newton = fp.read<double>();
        
        //----------------------------------------------------------------------
        // load fields
        //----------------------------------------------------------------------
        fp.load( h.entry, h.bytes);
        fp.load( w.entry, w.bytes);
        return true;
    }
    
    
    void loop(const size_t  num_iter, 
              const Real    dt, 
              const size_t  every,
              void (Simulation:: *process)(Real),
              const char *prefix)
    {
        eta ETA;
        const double eta_fac = 1.0 / num_iter;
        std::cerr.flush();
        assert(process);
        assert(prefix);
        initialize();
        
        if(GFX)
            gfxDraw();
        
        fprintf( stderr, "==== Starting Simulation type '%s'\n",prefix);
        const string perf_name = vformat("%s%u-perf.dat",prefix,unsigned(upper-1));
        save_perf( perf_name, 0, false );
        const string db_name   = vformat("%s%u.bin.gz",prefix,unsigned(upper-1));
        save_db( db_name, 0, false );
        
        double tmx_last = chrono.query();
        ETA.reset();
        for( size_t iter=1; iter <= num_iter; ++iter )
        {
            const Real t = dt * iter;
            (this->*process)(dt);
            if( 0 == (iter%every) )
            {
                const double tmx_ell = chrono.query() - tmx_last;
                const double fps     = every / tmx_ell;
                ETA(iter * eta_fac);
                const duration _done( ETA.time_done );
                const duration _left( ETA.time_left );
                fprintf( stderr, 
                        "  t = %10.3f @%7.0f fps in %02uh%02um%04.1f | ETA= %02uh%02um%04.1f     \r",
                        t,
                        fps,
                        _done.h,
                        _done.m,
                        _done.s,
                        _left.h,
                        _left.m,
                        _left.s); fflush( stderr );
                save_perf(perf_name, t, true);
                save_db(  db_name,   t, true);
                if(GFX)
                    gfxDraw();
                tmx_last = chrono.query();
            }
        }
        fprintf(stderr,"\n"); fflush(stderr);        
    }
    
    
    Real one_error( Real H, Real W ) const 
    {
        return Fabs(H*W-Kw)/Kw;
    }
    
    Real max_error() const throw()
    {
        Real ans = one_error(h[0],w[0]);
        for( size_t i=1; i <= N; ++i )
        {
            const Real tmp = one_error(h[i],w[i]);
            if( tmp > ans ) ans = tmp;
        }
        return ans;
    }
    
private:
    YOCTO_DISABLE_COPY_AND_ASSIGN(Simulation);
    void diffeq( array<Real> &dYdt, Real t, const array<Real> &Y )
    {
        const Real rate = (kd - k2 * Y[1] * Y[2]);
        dYdt[1] = rate; //! F[1]
        dYdt[2] = rate; //! F[2]
    }
    
    void jacobn( array<Real> &dFdt, matrix<Real> &dFdY, Real t, const array<Real> &Y )
    {
        dFdt[1] = 0;
        dFdt[2] = 0;
        
        dFdY[1][1] = dFdY[2][1] = -k2 * Y[2];
        dFdY[1][2] = dFdY[2][2] = -k2 * Y[1];
        
    }
    
};

#if defined(USE_FLTK)
#include <FL/gl.h>
#include <FL/Fl_Gl_Window.H>
#include <FL/Fl.H>

namespace {
    
    class MyWindow : public  Fl_Gl_Window
    {
    public:
        const Simulation &sim;
        MyWindow(int X, int Y, int W, int H, const char *L, const Simulation &s)
        : Fl_Gl_Window(X, Y, W, H, L),
        sim(s)
        {
            std::cerr << "** Creating OpenGL Window" << std::endl;
        }
        
        virtual ~MyWindow() throw()
        {
            std::cerr << "** Destroy OpenGL Window" << std::endl;
        }
        
        
        virtual void draw()
        {
            if (!valid()) 
            {
                valid(1);
                glLoadIdentity();
                glViewport(0,0,w(),h());
            }
            
            glClear(GL_COLOR_BUFFER_BIT);
            display( sim.X, sim.w, FL_GREEN  );
            display( sim.X, sim.h, FL_RED    );
        }
        
        inline void display( const Array1D &X, const Array1D &Y, Fl_Color c )
        {
            Real xmin = X[X.lower];
            Real xmax = X[X.upper];
            Real vmin=0,vmax=0;
            Y.get_min_max(vmin, NULL, vmax, NULL);
            Real ymin = -Log10(vmax);
            Real ymax = -Log10(vmin);
            const Real xlen = xmax-xmin;
            const Real dx   = 0.05 * xlen;
            const Real ylen = ymax - ymin;
            const Real dy   = 0.05 * ylen;
            glPushMatrix();
            glLoadIdentity();
            glOrtho(xmin-dx,xmax+dx, ymin-dy, ymax+dy, 1, -1);
            gl_color(c);
            glBegin(GL_LINE_STRIP);
            for( unit_t i=X.lower; i <= X.upper; ++i )
            {
                glVertex3d(X[i], -Log10(Y[i]), 0);
            }
            glEnd();
            
            glPopMatrix();
        }
        
        
    private:
        YOCTO_DISABLE_COPY_AND_ASSIGN(MyWindow);
    };
    
}

static  MyWindow *gfx = NULL;
static void gfxDraw()
{
    assert(gfx);
    gfx->redraw();
    Fl::check();
}
#else
static void gfxDraw() {}
#endif




static inline Real average_and_sigma( const array<Real> &v, Real &sig)
{
    Real sum = 0;
    const size_t n = v.size();
    for( size_t i=1; i <= n; ++i ) 
    {
        sum += v[i];
    }
    const Real ave = sum/n;
    sum = 0;
    for( size_t i=1; i <= n; ++i )
    {
        const Real tmp = v[i] - ave;
        sum += tmp * tmp;
    }
    sig = sqrt(sum/(n-1));
    return ave;
}



int main( int argc, char *argv[] )
{
    const char *progname = _vfs::get_base_name(argv[0]);
    try
    {
        //======================================================================
        // thread placements
        //======================================================================
        {
            const size_t np = hardware::nprocs();
            thread::assign_cpu( thread::get_current_handle(), np-1);
        }
        
        //======================================================================
        // activate graphical output or not
        //======================================================================
        {
            string GFX_VAL;
            const string NO_GFX = "NO_GFX";
            if( environment::get(GFX_VAL, NO_GFX) )
                GFX = false;
        }
        if( GFX )
            std::cerr << "-- GFX is activated" << std::endl;
        else 
            std::cerr << "-- GFX is disabled" << std::endl;
        
        //======================================================================
        //
        // Loading physical parameters from lua
        //
        //======================================================================
        if( argc <= 4)
            throw exception("usage: %s config.lua dt_save tmax tol", progname );
        
        //----------------------------------------------------------------------
        // create the lua virtual machine
        //----------------------------------------------------------------------
        Lua::State VM;
        L = VM();
        
        //----------------------------------------------------------------------
        // parse the file
        //----------------------------------------------------------------------
        {
            const string config_filename = argv[1];
            Lua::Config::DoFile(L,config_filename );
        }
        
        //======================================================================
        //
        // computing logical parameters
        //
        //======================================================================
        std::cerr << "-- Create Simulation" << std::endl;
        Simulation sim( L );
        Real dt_save = strconv::to_real<Real>( argv[2], "dt_save"); 
        Real tmax    = strconv::to_real<Real>( argv[3], "tmax" );
        sim.tol      = strconv::to_real<Real>( argv[4], "tol" );
        const string reaction = Lua::Config::Get<string>( L, "reaction");
        if( reaction != "explicit" && reaction != "implicit" )
            throw exception("invalid reaction '%s'", reaction.c_str());
        //======================================================================
        //
        // compute best dt
        //
        //======================================================================
        Real dt = 0.0;
        {
            const Real D_over_dx2 = max_of<Real>( sim.Dh, sim.Dw ) * sim.inv_dx * sim.inv_dx;
            const Real dt_max     = ALPHA_SAFE / D_over_dx2;
            std::cerr << "-- dt_max = " << dt_max << std::endl;
            const Real dt_log      = pow(10.0,floor(log10(dt_max)));
            dt = floor(dt_max/dt_log) * dt_log;
            std::cerr << "-- dt=" << dt << std::endl;
        }
        
        //======================================================================
        //
        // compute iteration parameters
        //
        //======================================================================
        if( tmax < dt ) tmax = dt;
        size_t       num_iter = Ceil(tmax/dt);
        const size_t every    = (num_iter * dt_save) / tmax;
        while( 0 != (num_iter%every) ) ++num_iter;
        std::cerr << "-- tmax      = " << tmax     << std::endl;
        std::cerr << "-- dt_save   = " << dt_save  << std::endl;
        std::cerr << "-- num_iter  = " << num_iter << std::endl;
        std::cerr << "-- every     = " << every    << std::endl;
        std::cerr << "-- Tolerance = " << sim.tol  << std::endl;
        
        //======================================================================
        //
        // Loop
        //
        //======================================================================
#if defined(USE_FLTK)
        auto_ptr<MyWindow> win;
        if(GFX)
        {
            const int MySize = 512;
            win.reset( new MyWindow(10,10,MySize,MySize,"Simulation",sim) );
            gfx = win.__get();
            win->mode( FL_RGB );
            win->show();
            Fl::check();
        }
#endif
        sim.loop( num_iter, dt, every, & Simulation:: relaxed_scheme        , "ch");
        sim.loop( num_iter, dt, every, & Simulation:: implicit_plain_scheme , "im");
        
        
        //======================================================================
        //
        // database processing
        //
        //======================================================================
        std::cerr << "-- Processing Databases" << std::endl;
        {
            Simulation &sim_ch = sim;
            Simulation  sim_im( L );
            const unsigned nv = sim_ch.upper-1;
            
            ios::igzstream fp_ch( vformat("ch%u.bin.gz", nv) );
            ios::igzstream fp_im( vformat("im%u.bin.gz", nv) );
            
            Real t=0;
            Real ch_max_error = 0;
            Real im_max_error = 0;
            Real im_max_norm  = 0;
            vector<Real> im_speedup;
            const size_t num  = sim_ch.X.width;
            vector<Real> X(  num,as_capacity );
            
            vector< vector<Real> > pH(32,as_capacity);
            vector<Real>           pH_time(32,as_capacity);
            const Real             frame_time = 60.0;
            Real                   last_frame = -2*frame_time;
            
            for( unit_t i=sim_ch.X.lower;i<=sim_ch.X.upper;++i)
                X.push_back( sim_ch.X[i] );
            
            {
                ios::ocstream  fp( vformat("results%u.dat",nv), false );
                fp("#t err_ch err_im norm_im t_diff_ch t_diff_im t_proj_ch t_proj_im speedup ave_newton\n");
                unsigned frame_id = 0;
                while( sim_ch.load_db(t, fp_ch ) )
                {
                    fprintf( stderr, "Process Frame %6d @t= %10.4f        \r", frame_id, t); fflush( stderr );
                    Real tmp = 0;
                    if( !sim_im.load_db(tmp,fp_im) )
                        throw exception("Missing implicit frame");
                    
                    //----------------------------------------------------------
                    // store this frame ?
                    //----------------------------------------------------------
                    if( t - last_frame >= frame_time )
                    {
                        pH_time.push_back(t);
                        {
                            const vector<Real> tmp; pH.push_back(tmp);
                        }
                        vector<Real> &tgt = pH.back();
                        for( unit_t i=sim_ch.X.lower;i<=sim_ch.X.upper;++i)
                            tgt.push_back( -Log10(sim_ch.h[i]) );
                        last_frame = t;
                    }
                    
                    
                    
                    const Real ch_total       = sim_ch.t_diff+sim_ch.t_proj;
                    const Real im_speedup_val = ch_total > 0 ? (sim_im.t_diff+sim_im.t_proj)/ch_total : 1;
                    
                    if( ch_total > 0 )
                    {
                        im_speedup.push_back(im_speedup_val);
                    }
                    
                    const Real ch_max_error_val = sim_ch.max_error();
                    const Real im_max_error_val = sim_im.max_error();
                    
                    if( ch_max_error_val > ch_max_error ) ch_max_error = ch_max_error_val;
                    if( im_max_error_val > im_max_error ) im_max_error = im_max_error_val;
                    
                    Real im_norm = 0;
                    for( unit_t i=sim_ch.X.lower;i<=sim_ch.X.upper;++i)
                    {
                        const Real h_ch = sim_ch.h[i];
                        const Real w_ch = sim_ch.w[i];
                        const Real h_im = sim_im.h[i];
                        const Real w_im = sim_im.w[i];
                        const Real dh   = h_ch - h_im;
                        const Real dw   = w_ch - w_im;
                        const Real tmp  = max_of<Real>( Fabs(dh)/Fabs(h_ch+h_im), Fabs(dw)/Fabs(w_ch+w_im) );
                        if( tmp > im_norm ) im_norm = tmp;
                    }
                    
                    if( im_norm > im_max_norm ) im_max_norm = im_norm;
                    
                    fp( "%g",  t);
                    fp( " %g", ch_max_error_val );
                    fp( " %g", im_max_error_val );
                    fp( " %g", im_norm      );
                    fp( " %g", sim_ch.t_diff);
                    fp( " %g", sim_im.t_diff);
                    fp( " %g", sim_ch.t_proj);
                    fp( " %g", sim_im.t_proj);
                    fp( " %g", im_speedup_val);
                    fp( " %g", sim_ch.ave_newton);
                    fp("\n");
                    ++frame_id;
                }
                fprintf( stderr, "\n");
                
            }
            
            
            std::cerr << "#pH =" << pH.size() << "/" << pH_time.size() << std::endl;
            std::cerr << "pH@" << pH_time << std::endl;
            {
                ios::ocstream fp( vformat("profile%u.dat",nv), false );
                fp("#t");
                for( size_t i=1; i <= pH_time.size(); ++i )
                {
                    fp(" ph@%g", pH_time[i]);
                }
                fp("\n");
                for( unit_t i=sim_ch.X.lower;i<=sim_ch.X.upper;++i)
                {
                    const size_t j = 1 + (i-sim_ch.X.lower);
                    fp("%g", sim_ch.X[i] );
                    for( size_t k=1; k <= pH.size(); ++k )
                    {
                        fp(" %g", pH[k][j]);
                    }
                    fp("\n");
                }
            }
            
            {
                ios::ocstream fp( vformat("info%u.dat",nv), false );
                Real       speedup_sig_im = 0;
                const Real speedup_ave_im = average_and_sigma( im_speedup, speedup_sig_im);
                const Real total_time_ch  = sim_ch.t_diff + sim_ch.t_proj;
                const Real total_time_im  = sim_im.t_diff + sim_im.t_proj;
                fp("%g %g %g %g %g %g %g\n",
                   Log10(sim.tol), 
                   Log10(ch_max_error), 
                   Log10(im_max_error),
                   total_time_ch,
                   total_time_im,
                   speedup_ave_im, 
                   speedup_sig_im);
            }
        }
        
        
    }
    catch( const exception &e )
    {
        std::cerr << "******** " << e.what() << std::endl;
        std::cerr << "******** " << e.when() << std::endl;
        return -1;
    }
    catch(...)
    {
        std::cerr << "******** Unhandled exception!" << std::endl;
        return -1;
    }
    return 0;
}
